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Abstract. - We propose a simple phenomenological model for wet granular media to take into 
account many particle interaction through liquid in the funicular state as well as two-body cohesive 
force by a liquid bridge in the pendular state. In the wet granular media with small liquid content, 
liquid forms a bridge at each contact point, which induces two-body cohesive force due to the 
surface tension. As the liquid content increases, some liquid bridges merge, and more than two 
grains interact through a single liquid cluster. In our model, the cohesive force acts between 
the grains connected by a liquid-gas interface. As the liquid content increases, the number of 
grains that interact through the liquid increases, but the liquid-gas interface may decrease when 
liquid clusters are formed. Due to this competition, our model shows that the shear stress has a 
maximum as a function of the liquid-content. 



Introduction. — It is well known that just adding 
small amount of liquid to dry granular materials changes 
their behaviors drastically [1~3]. Liquid forms bridges be- 
tween grains, and the bridge induces cohesion between 
grains due to the surface tension. When the amount of liq- 
uid is small, liquid bridges are formed between two grains 
at contact points; such a state is called the pendular state 
(Fig. [Ha)). As the amount of liquid is increased, several 
liquid bridges are connected to form a liquid cluster that 
, includes more than two grains as in Fig. \V[h) , which is 
' called the funicular state. As we add the liquid further, 
all the pores will be filled with the liquid. It is called 
the capillary state if the liquid pressure is still lower than 
the air pressure, and finally the system becomes the slurry 
state if the amount of liquid is enough to make the liquid 
pressure as high as the air pressure. 

The liquid content dependence of the mechanical prop- 
erties of wet granular materials has not been well under- 
stood yet. The most studied case is the pendular state, 
both experimentally and numerically [1-6,6-9]. There, 
the situation is rather simple because the grain-grain in- 
teractions, including the cohesion due to liquid bridges, 
are mainly two-body interactions. However, as the liquid 
content is increased, liquid clusters become non-negligible, 
which induce many-body interaction via the liquid. Meth- 
ods and/or models to deal with such a situation have 
not been established yet and still under active research 




Fig. 1: Schematic description of the distribution of liquids in 
the wet granular materials, (a) The pendular state, where 
liquid bridges are formed at contact points, (b) The funicular 
state, where some liquid clusters that connect more than two 
grains are formed. 



( [10-15]). 

Recent experiments show that, upon increasing the liq- 
uid content, there is a maximum in the material yield 
stress at rather low liquid content, in the pendular state or 
the funicular state [10-12]. Especially, M0ller et al. [12] 
showed a clear increase of the shear modulus with the 
liquid content in the pendular state and a drop in the fu- 
nicular state. 

In this paper, we propose a simple phenomenological 
model that takes into account the many-body interaction 
through the liquid clusters in the funicular state. In our 
model, the cohesive force acts among the grains connected 
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by the liquid-gas interface. As the liquid content increases, 
the number of grains that interact through the liquid in- 
creases, but the liquid-gas interface may decrease when the 
liquid cluster is formed. We show that this competition 
results in a maximum in the shear stress as a function of 
liquid-content, which may correspond qualitatively with 
the experimentally observed behavior of the shear modu- 
lus as a function of liquid content. 

Model. — In this section, we first introduce a model 
for the pendular state, where the two-body cohesion due 
to liquid bridges is incorporated in the soft sphere model 
of dry grains with linear elastic force and dissipation (Dis- 
crete Element Method, DEM [16]). For the cohesive inter- 
action, history-dependent formation of liquid bridge pro- 
posed by Schultz et al. [17] is adopted. We then extend the 
model to the funicular state, where some pores are filled 
with liquid and the many-body interactions acts among 
particles. 

Model P for the pendular state with two-body inter- 
action. For simplicity, we consider a frictionless two- 
dimensional system, where grains are modeled by poly- 
disperse disks. Let us consider the interaction between 
particles i and j of diameters ai and (jj and masses 
and TTij, at positions and rj with velocities Vi and Vj, 
respectively. There is no interaction between them before 
they collide. Once they are in contact, i.e., |ry| = \ri — rj\ 
becomes less than (di aj)/2, the force on the particle i 
by j is given as a function of A = (ui -f (Jj)/2 — \rij \ by 

F,, = /(A)n,, - 7j{A)v,, (1) 

with the normal unit vector riy ~ f'ij/l'f'ijl the rela- 
tive velocity Vij = Vi — Vj, as long as \rij \ < a{ai -{-a-j)/2; 
The parameter a{> 1) determines the distance between 
grains at which the liquid bridge breaks, and once \rij\ 
exceeds a{ai -\- aj)/2, Fij = and the pair i,j do not 
interact until they are in contact again. 

The functions /(A) and r7(A) in the interaction of 
Eq.([T]) are given by 

/(A) = max(fcA,-i^o) (2) 

and 

^(A) = I ^s-in for A > 0, 

' \ ?yiiquid for A < 0. ^ ' 

We can interpret that A characterizes the overlap (for 
A > 0) or the spacing (for A < 0) between the grains. 
The parameter k characterizes the elastic constant, 77(A) 
is the damping parameter, and a positive constant _Fo char- 
acterizes the cohesion due to a liquid bridge 0. 

^ One might presume that the cohesive force should act even 
when grains are in contact. This can be dealt with, for the present 
frictionless system, simply by reinterpreting A as being measured 
from the overlap position where the cohesive force and repulsive 
force due to grain elasticity balance. The dissipation term of eq.(3) 
should be modified accordingly, but such a minute adjustment barely 
affects the overall system behavior. 
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Fig. 2: Schematic description of interaction of (a) the model 
for the pendular state (Model P) and (b) the model for the fu- 
nicular state (Model F). The solid circles show grains and the 
dotted line indicates the interaction range after the formation 
of a liquid bridge. The solid arrows indicates the pair of grains 
that cohesive force acts. The dotted arrow in (b) shows "inac- 
tive bonds", i.e., they are interacting via liquid cluster but the 
bond exert no cohesive force. 

The damping parameter 77(A) in eq.([3]) depends on A 
because the dissipation when grains are in contact (char- 
acterized by 77grain) and the dissipation via liquid viscosity 
(characterized by jyiiquid) has different physical origin and 
their values can be quite different. 

Notice that the model with a — \ corresponds to dry 
granular materials with neither cohesive force nor hystere- 
sis, and a is greater than unity for wet cohesive grains. 
When a > 1, there is a hysteresis in the interaction [17], 
and this hysteresis causes dissipation in addition to that 
caused by the viscous force. 

In this model for the pendular state (here we call it as 
Model P), all interactions are given as two-body forces. 
The situation is schematically illustrated in Fig. [2Ua); the 
attractive force between the particles connected by solid 
arrows represents the cohesion due to liquid bridges. 

Model F for the funicular state with many-body inter- 
action. When the amount of liquid is increased enough 
for the system to be in the funicular state, some pores are 
filled by liquid, and more than two grains are connected by 
a single liquid-filled region. Let us consider the situation 
shown in Fig.[2jb), where four particles are connected by a 
liquid cluster. In this case, there is no liquid-gas interface 
that connects particles 2 and 4; in contrast to the pendular 
state (Fig. [2ja)), there is no direct cohesive force between 
them, though the cohesive force between the particles con- 
nected by solid arrows in Fig. Eljb) still tend to hold all 
the particles together. There should be still some effec- 
tive cohesive force between the particles 2 and 4 caused 
by the lower pressure in the liquid phase than that in the 
air phase, but here we disregard this. 

We model such a multi-particle effect by making the 
bond inactive when the bond connects the two particles 
between which a liquid cluster exists. The presence of the 
liquid cluster between two grains is determined by the sim- 
ple criterion: The grains i and j are connected by a liquid 
cluster without liquid-gas interface if they are connected by 
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Fig. 3: Snapshot of grain-grain interactions. The interacting 
grains are connected by bonds. The cohesive bonds, or the 
bonds with Ni^j < nth and A < 0, are shown by red solid lines, 
the repulsive bonds with liquid interfaces, or Nij < nth and 
A > 0, are by blue long-dashed lines, the repulsive bonds in 
liquid clusters, or Nij > nth and A > 0, are by green long- 
dashed lines, and the inactive bonds, or Ni,j > nth and A < 0, 
by a light blue short-dashed line. 



a liquid bridge and the number of grains Nij that are con- 
nected to both i and j by liquid bridges is larger or equal 
to a certain threshold value nth- The liquid bridge is the 
same as defined in the previous subsection, and we take 
nth = 2 for the two-dimensional system. If two grains are 
connected by a big cluster without liquid-gas interface, 
the liquid bridge between them is inactive and no cohe- 
sive force acts through. An example in two-dimensional 
system with nth = 2 is given in Fig. [Hb) . 

With this criterion, the force on grain i by j is given not 
by Eq. ^ but by 



/(A) [e(/(A)) 

-??(A)t»y, 



(1 - e(/(A)))eK;, - iv,,,)] 



(4) 



where /(A) is given by Eq. ([2]) and Q{x) is the Heav- 
iside step function defined as Q{x) = 1 for a; > and 
@{x) =0 for X < 0. We call this model for the funicular 
state as Model F here. Note that it is straightforward to 
extend this model to a three-dimensional system by choos- 
ing nt/i = 3. 

Figure [3] illustrates an example of the interaction net- 
work in Model F. The interacting grains are connected by 
bonds. The cohesive bonds, or the bonds with Nij < nth 
and A < 0, arc shown by red solid lines, the repulsive 
bonds with liquid interfaces, or Nij < nth and A > 0, 
are by blue long-dashcd lines, the repulsive bonds in liq- 
uid clusters, or Nij > nth and A > 0, are by green 
long-dashed lines, and the inactive bonds, or Nij > nth 
and A < 0, by a light blue short-dashed line. We can 



Fig. 4: Schematic description of the simulated configuration. 



see several clusters of particles connected by bonds with 
> nth, which corresponds to liquid clusters. 

In this model, the parameter a for the interaction range 
represents the liquid content. In Model P, increasing a 
just increases the number of cohesive bonds. In Model 
F, on the other hand, not only the number of interacting 
particles increases, but also the fraction of bonds that are 
surrounded by other bonds {Nij > nth) increases; such 
bonds exert no cohesion. Thus, the larger a corresponds 
to the case with more liquid content. 

We now examine how the competition between these 
effects of increasing a in Model F affects the response of 
the wet granular assembly against shear. 

Simulation Setup. — Wc mainly focus on Model F 
under shear, but we also present the results of Model P 
for comparison. 

The grains are modeled by polydisperse disks whose di- 
ameters are uniformly distributed between 0.8cr to cr, and 
all the grains have the same mass m. The system is peri- 
odic in the x direction. In the y direction, there are two 
parallel rough walls of mass M, each of which consists of 
N.j^ particles with diameter cr glued without spacing (Fig. 
m the total length of the wall L = aN^^). 

The top (bottom) wall is moving in the x direction at 
a velocity U {—U) in order to realize a shear flow with 
velocity gradient in the y direction (Fig. S}. A constant 
pressure is applied to the walls in the y direction, and 
the y coordinate of the top (bottom) wall Yt (!&) obeys 
the equation of motion MYt = Fw,t{t) — PwL (MYi, = 
Fw,b{t) + PwL), where F^,f (i) (i^u,,6(t)) is the y-component 
of the force exerted from particles on the upper (bottom) 
wall at time t. 

For Model F, we also need to decide whether there is 
liquid bridges among grains that consist the wall; That 
affects the formation of the liquid cluster near the wall. 
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We initially impose that a neighboring pair of grains on 
the wall has a liquid bridge with the probability 1/2, and 
the bridges among wall grains do not disappear or newly 
created during the simulation. 

Initial configurations for simulations are prepared as 
follows; N particles are distributed randomly with large 
enough spacing, and the wall pressure is applied with- 
out shear {U — 0). After waiting long enough so that all 
the grains are at rest, we start to shear the assembly at a 
constant U. 

The data are collected in the steady states, where 
both the average kinetic energy of grains and the the y- 
coordinates of the walls become almost constant. 

The parameters are given in the dimensionless form with 
the unit length a, imit mass m, and the unit time scale 
y/ lOOma'^ / Fq in the following. In this unit, we choose the 
elastic constant fc to be 5 x 10'*, the viscous coefficient due 
to grain contact lygrain to be 10 (which gives the restitution 
coefficient e„ for dry grains (a—l) around 0.9), and the 
viscous coefRcient due to liquid bridge Tyuquid to be 1. In 
this paper, we focus on the situation where the relative 
velocity of particles is 0(1), thus the viscous force is about 
1/100 of the cohesive force. Also, note that the particle 
deformation is typically less than 1% of its diameter under 
the force applied in the simulations. 

We use the wall mass M — N^^m, and several values of 
the wall pressure Pw The system size are mostly L — AO 
and the number of particles in the system N — 1600 un- 
less otherwise noted. In most of the cases, the wall ve- 
locity U is chosen to give the shear rate 7 « 1 in the 
middle of the system, and smaller/larger U are also exam- 
ined to see the tendency. The liquid-content dependence is 
studied by changing the interaction range a. The second- 
order Adams-Bashforth method and the trapezoidal rule 
are used to integrate the equations for the velocity and 
position, respectively, with a time step for integration 
dt = 10-4. 

Simulation Results. — We will confirm that the sys- 
tem can achieve plain shear flows first, and then we will 
investigate the effect of liquid content on the model be- 
havior by changing the interaction range a. 

First we check whether a plain shear flow can be realized 
in the present simulation. Figure [5] (a) and (b) shows the 
velocity profile and the packing fraction for Model F with 
U = 18, Pu, = 100, and a = 1.3. From Fig. EKa), we see 
that the realized shear rate is almost constant with 7 « 1, 
whereas some spatial structure can be seen in the packing 
fraction in Fig. EJb). Similar profiles has been obtained 
with various a and Pu, > 50 for Model F. However, we 
found that, when the wall pressure P^ is too weak and/or 
the system size is too big, most of the grain may stick to 
one of the walls because of the cohesion, and the velocity 
gradient is confined close to the other wall. This tendency 
is stronger for Model P with large a and large enough 
system (FiglSJc) solid line (L=40)). When the system is 
small enough (Fig. [5Kc) dashed line (L=20)), the plain 
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Fig. 5: (a) and (b): The velocity (a) and the packing fraction 

(b) profile of the system with U — 18 for Model F with a — 1.3 
and Pw = 100. In (a), we can see that the shear rate 7 is almost 
constant over the system with 7 « 1, while a spatial structure 
can be seen in the packing fraction profile near the boundary. 

(c) and (d): The velocity (c) and the packing fraction (d) profile 
of the system for Model P with a = 1.3 and P„ = 100. The 
data for L = 40 and [/ = 18 is shown by solid lines, and the 
data for L = 20 and (7 = 9 is shown by dashed lines. 



shear can be obtained for Model P also. In the following 
we focus on the parameter range where plain shear flow 
can be realized approximately. 

The interaction range a dependence of the shear stress 
S is shown in Figure [6] in the case of P^ = 100 for Model 
F and Model P with 7 « 1. (for Model F, U = 18 is 
used, while for Model P the small system with L — 20 
and N = 400 with [/ = 9 is studied to avoid the non- 
uniform shear). Note that, for a = 1 (no cohesion due to 
liquid bridge), both models should give the same results. 
In Model P, the shear stress increases with a, because 
larger a simply results in more cohesive interactions. In 
Model F, on the other hand, S shows a maximum around 
a = 1.2. 

The maximum of the shear stress S* as a function of a in 
Model F can also be seen for various values of pressure P^ , 
but the location of the maximum changes as can be seen 
in Fig. [ZKa). The packing fraction cj) at the middle of the 
system, is plotted against a for various P^ in Fig. [7]^b); 
The packing fraction is higher for larg expected, 
but the dependence on a is rather weak. 

The existence of the maximum in the shear stress should 
be due to the competition of the increasing interacting par- 
ticles and the liquid-cluster formation which reduces the 
cohesive interaction. To see this more directly, the aver- 
age number of interacting particles ( "bonds" ) per particle 
Nint and the number of cohesive bonds (i.e., the bond 
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Fig. 6: The interaction range a dependence of the shear stress 
S in the case of P„ = 100 for Model F and Model P with 
7 « 1. (for Model F, [/ = 18 is used, while for Model P the 
small system with L — 20 and A*' = 400 with U — 9 is studied 
to avoid the non- uniform shear). For a = 1, both models give 
the same results. 

that satisfies Nij < nth and A < 0) per particle Nc are 
shown in Fig.[71Jc) and (d), respectively. We can see that, 
even though the interacting particles increase with a, the 
cohesive bonds has a maximum as a function of a. The 
location of the maxima are found to be slightly different 
between the shear stress S and the number of cohesive 
interactions but they have a similar trend that the 
maximum moves towards smaller a for larger P^, . We also 
found that is actually larger for smaller P^, although 
Nint is larger for larger P„,. This can be understood as 
an effect of the packing fraction </), i.e., for larger P^, the 
packing fraction cf) increases (Fig. [7}d), which results in 
larger Nint but Nc decreases because many of the cohe- 
sive interactions become inactive as Nj tends to exceed 
nth more easily. 

The decrease of the cohesive bond upon increasing a 
should be a direct result of the increasing size of liq- 
uid clusters. We characterize the liquid cluster size by 
Ncius, the number of particles connected by bonds with 
Nij > nth- The cluster size distribution P{Ncius) is shown 
in Fig. [5] for = 100 and 7 w 1. We find that the distri- 
bution show exponential decay for small a, but the decay 
becomes slower as a increase, and shows power-low like 
behavior at around a w 1.4, beyond which we get a finite 
probability to find a cluster as big as the system size. The 
cluster size distribution depends on the pressure and 
the shear rate 7, but the tendency to have an exponential 
decay for small a and the slower decay for larger a is seen 
for all the cases (data not shown). 

Summary and Discussion. — We proposed a simple 
model for the wet granular media in the funicular state. 
We assumed that the cohesive force acts between the par- 
ticles connected by the liquid-gas interfaces but not those 
in a liquid cluster. We examined the model behavior upon 
changing the liquid content, which is represented by the 




Fig. 7: The interaction range a dependence of the shear stress 
S (a), the packing fraction (j) (b), the number of bonds per 
particle Ni„t (c), and the number of cohesive bonds per particle 
Nc (d) in Model F for various pressure P^. U — 18, which gives 
7«1. 

interaction range parameter a. The shear stress S* in a 
steady shear flow was demonstrated to show a maximum 
at a certain value of a. The existence of the shear stress 
maximum is a result of competition between the two ef- 
fects upon increasing the liquid content, i.e., the increase 
of interacting bonds, which increases the shear stress, and 
the formation of liquid clusters, which reduces the cohesive 
interaction, thus reduces the shear stress. 

In experiments, the shear modulus has been demon- 
strated to show a maximum upon increasing the liquid 
content [12]. The measurements have been done for small 
enough stain at finite frequencies using rheometer. The 
observed shear modulus is fairly independent of the fre- 
quency in the range between 0.01 Hz and 10 Hz, and found 
to show maximum values around the transition regime be- 
tween the pendular and the funicular regimes for various 
granular materials. This experimentally observed maxi- 
mum in the shear modulus may be interpreted as a result 
of the competition between the formation of liquid bonds 
and liquid clusters discussed above. 

Another example of experimental observation is on the 
rotating drum [9] . Within the pendular state, it has been 
shown that the surface angle increases as the liquid vol- 
ume increases at the low rotation rate, while at the high 
rotation rate the angle decreases with the liquid volume. 
The behavior at the low rotation rate is consistent with 
that of our Model P, but that at the high rotation rate is 
not. The latter case seems to indicate that the lubrication 
and viscous forces are important at higher rotation rate, 
and the rheological properties may not be simulated by 
simple two-body cohesions. 
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Fig. 8: The cluster size distribution P{Ncius) with Pw = 100 
and 7 Ri 1 for various interaction range a. 

Before concluding, let us discuss limitations and possi- 
ble extensions of our model. 

Firstly, the lubrication effect is not included in the 
present model, and that will be important for high shear 
rate cases as mentioned above. 

Secondly, the present model should not be valid for 
both limiting cases of the small liquid content and the 
fully saturated regime. In the former case, the friction 
should become important for both the sliding [5] and the 
rolling [18-20] motions because the cohesive force pulls 
the particles together to contact even without external 
pressure while the liquid lubrication suppresses the fric- 
tion effect for larger liquid content. In the fully saturated 
regime, the rheology will be determined by the hydrody- 
namic interaction, which is not included in the present 
model. 

Lastly, the liquid content is only represented by the 
parameter a, and its conservation is not taken into ac- 
count even though this may affect the configuration of the 
interaction network. The model for the pendular state 
proposed by Richefeu et al. [7,8] is one of the possible ap- 
proach to include the liquid conservation effectively, where 
a liquid volume dependent interaction has been adopted 
and the liquid is redistributed to each liquid bridge ac- 
cording to a simple rule along the time evolution. This 
model showed the saturation of cohesion effects as a func- 
tion of the liquid content. Thus, if one add the liquid 
cluster criterion to it as in our model, the drop of the co- 
hesion effect at high liquid content will be also observed 
because the liquid cluster will reduce the cohesive bonds. 
It will be interesting to see how the behavior of the model 
changes by including these effects. 
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